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In light of recent proposals to realize a topological superconductor on the surface of strong topo¬ 
logical insulators, we study impurity and vortex scattering in two dimensional topological super¬ 
conductivity. We develop a theory of quasiparticle interference in a model of the surface of a three 
dimensional strong topological insulator with a pairing term added. We consider a variety of dif¬ 
ferent scatterers, including magnetic and nonmagnetic impurity as well as a local pairing order 
parameter suppression associated with the presence of a vortex core. Similar to the case of a surface 
of a three dimensional topological insulator without pairing, our results for non-magnetic impurity 
can be explained by the absence of back scattering, as expected for a Dirac cone structure. In 
the superconducting case, doping away from the Dirac point leads to a doubling of the contours 
of constant energy. This is in contrast to the unpaired case where the chemical potential simply 
adds to the bias voltage and shifts the energy. This doubling of contours results in multiplying the 
number of possible scattering processes in each energy. Interestingly, we hnd that some processes 
are dominant in the impurity case while others are dominant in the vortex case. Moreover, the two 
types of processes lead to a different dependence on the chemical potential. 


I. INTRODUCTION 


The search for a Majorana mode in a condensed mat¬ 
ter system is currently a very active area of research. 
Some very promising experimental results have been re¬ 
ported in one dimensiorP^ while in two dimensions ap¬ 
pealing theoretical proposals exislP^. The above sys¬ 
tems are based on the combination of momentum spin 
locking with superconducting pairing. In two dimensions 
this combination is provided by heterostructures of a su¬ 
perconducting layer in contact with either a three di¬ 
mensional topological insulator (3DSTI)P or a spin-orbit 
coupled semiconducto]®^. While the superconductor fur¬ 
nishes the pairing, the spin orbit coupled layer provides 
the momentum-spin lockin^^HUl^ Another type of pro¬ 
posal makes use of an innate tendency for developing 
superconductivity in materials with spin orbit coupling 

(socpuni. 

Given the above, it is timely to search for unique prop¬ 
erties of a system where superconductivity arises from an 
underlying Dirac-like band structure. In this Paper we 
develop a theory of quasiparticle interference (QPI) pat¬ 
terns in such systems. The QPI patterns are of direct 
experimental relevance as they can be measured using 
the technique of Fourier-transform scanning tunnelling 
spectroscopy (FT-STS)P2®^. This method measures the 
local density of states (LDOS) of a sample in the vicin¬ 
ity of a single impurity. Theoretically, the pattern which 
is observed is dictated by the underlying clean system 
{i.e without the impurity) and so properties of the clean 
system can be deduced from such a pattern and readily 
compared with calculation^^^H^. 

We consider the surface of a strong topological insu¬ 
lator with pairing added. We do this by treating two 
complementary models. The first is a single Dirac cone 
in the continuum and should capture universal proper¬ 


ties of a Dirac band structure in the presence of pairing. 
Second, we scrutinize our continuum model results by 
looking at a more physical lattice model. This model has 
been developecP^ for the surface of a strong topological 
insulator. It specifically avoids the doubling theorem by 
including both surfaces of the material. We calculate the 
QPI pattern for a variety of scatterers, including a charge 
defect, a magnetic defect and a defect in the supercon¬ 
ducting order parameter. 

The QPI for a strong topological insulator in the pres¬ 
ence of charge and magnetic impurities has been calcu¬ 
lated by one of us in the pasfP^. With the chemical po¬ 
tential tuned to the Dirac point, we obtain results con¬ 
sistent with this previous study with the exception that 
QPI pattern is only observable for energies above the su¬ 
perconducting gap Aq, as illustrated in Fig. |Ij As noted 
previously, unlike in normal metals where QPI shows dis¬ 
tinctive peaks, patterns in STIs with non-magnetic im¬ 
purities are non-singular and exhibit only an edge. We 



FIG. I: Quasiparticle interference patters, observable by FT- 
STS, on the surface of a strong topological insulator with 
SC order. Panels a) and b) show the effect of non-magnetic 
impurities on the normal (Ao = 0) and SC state, respectively. 
Panel c) displays the effect of disorder in the order parameter 
amplitude. 
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find this behavior to persist when the SC order is in¬ 
cluded. An interesting difference however occurs when 
we consider disorder in the SC order parameter ampli¬ 
tude A which will generically arise in the presence or 
pair-breaking impurities of vortices. In this case, QPIs 
display a sharp peak, as illustrated in Fig. [I]:). 

Moving the chemical potential away from the Dirac 
point does not result in any angular dependence of the 
properties of this pattern. However, we do see some in¬ 
teresting behavior in these patterns. We find a contour 
of singular LDOS features. As the chemical potential is 
increased (decreased) with respect to the Dirac point, the 
radius of the contour in the impurity induced QPI pat¬ 
terns is also increased (decreased) and changes linearly 
with the chemical potential. On the other hand, in the 
order parameter suppression QPI pattern, the radius of 
the singular contour is independent of the chemical po¬ 
tential. We trace this difference back to the fact that 
different quasiparticle scattering processes are favored in 
these two situations. It can also be argued that this par¬ 
ticular dependence on fi for these two types of scatterers 
comes from the underlying Dirac band structure. These 
observations in the continuum model can be understood 
heuristically via chirality angle arguments and are sup¬ 
ported by our lattice model calculations. 

Although tuning the chemical potential in these sys¬ 
tems may be challenging from an experimental point of 
view, the underlying result {e.g. the behavior of the QPI 
singularity) can be observed via other means. As an ex¬ 
ample, we suggest a simplified scenario where the sample 
we consider is placed in a capacitor. This adds a bias volt¬ 
age between the edges of the material. The voltage acts 
as a chemical potential with differing signs on each edge 
of the sample. By changing the bias one can observe the 
distinct behavior the QPI due to an impurity or a vortex. 
There has recently been an experimental realization of a 
set-up analogous to this capacitor systenP^. In this work 
it was shown that the chemical potential on each edge of 
thin film 3DSTI can be tuned independently of the other 
edge through the use of dual-gate structures. Our results 
are then directly applicable to an FT-STS measurements 
on such a device when the chemical potential on one edge 
is tuned to be opposite to each other edge. 

The rest of this paper is organized as follows. In the 
next section we give an overview of the methods we have 
followed to obtain our results. We present our method 
for calculating the QPI pattern, define each type of per¬ 
turbation, present our two model systems and discuss 
numerical details. Section m then presents our results 
for the QPI patterns. We begin by tuning the chemi¬ 
cal potential to the Dirac point, where the QPI pattern 
for our continuum model can be calculated exactly (the 
details of this are relegated to the Appendix). This not 
only helps us compare our results with a previous study, 
but also allows us to develop some intuition. We then 
move on to the main result of this paper, the behavior 
of QPI patterns at finite chemical potential. Here we 
present numerical data as well as further discussion. We 


then discuss an alternative to tuning the chemical poten¬ 
tial and perform an explicit calculation using our lattice 
model. 


II. METHODOLOGY 


A. Green’s function for a single impurity and the 
Born approximation 


We begin by discussing a general expression for the 
total Green’s function. The approach taken here will be 
similar to previous worlP^, but appropriately generalized 
to include superconductivity and to be suitable for ap¬ 
plication to our lattice model. In either case we take a 
total Hamiltonian H = Hq + Himp, the first term de¬ 
scribes the underlying clean system while the latter term 
describes the perturbation. The clean system’s Hamilto¬ 
nian is given by: 

= ( 1 ) 

k 


here we have defined = (ck,a, where a labels 

additional degrees of freedom in the model (spin, sub 
lattice, etc.). The matrix ?^k is given by 

/Mk) A(k) \ 

VAt(k) -ht(-k); 

where Ak and h(k) are matrices in the space of the de¬ 
grees of freedom implicit in a before, h(k) describes the 
normal band structure of the model and A(k) describes 
the pairing. Next we consider a local perturbation which 
couples to the electronic degrees of freedom in some gen¬ 
eral form: 


Himp = ^ E - kOV’k' (3) 

^ k,k' 

where N is the number of lattice sites; i/^^^(k — k') is 
the Fourier transform of the perturbation potential. Its 
matrix structure depends on the impurity we consider 
and will be discussed in detail below. 

Next we consider calculating the Green’s function for 
this model. We define the Green’s function as follows 


G(k,k',r) = -(T, 


^k(r)t/>j,,(0) ) 


(4) 


where r is the imaginary time and '0k(^) is the Heisen¬ 
berg picture version of ^/^k- Fourier transforming to Mat- 
subara frequency and defining the T-matrix we write 

G{k,k',iujn) = G°{k,iujn)Sk,k' (5) 

+ G°{k,iuJn)T{iuJm,k,k')G°{k',iuJn) 

where G°(k, zw„) = — is the Green’s func- 

tion of the clean system. In the case of a local perturba¬ 
tion, where the potential is a (^-function in real space, the 
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potential and the T-matrix are momentum independent 
and the self consistency equation for the T-matrix can 
be solved: 

T{iw^) = AT-iylp - AT-i ^ G°(k, icon)V°^p^ . 

(6) 

In principle poles in the T-matrix reveal information 
about bound states in the system. However, in this 
study we are interested in the QPI patterns in momen¬ 
tum space. Since the T-matrix has no momentum depen¬ 
dence, the results of a weak perturbations will not differ 
from the simpler Born approximation. To this end, we 
keep only the leading order term in which leaves 

G{k,k',iuJn) ^ G0(k,zw„)(5k,k' (7) 

+ - k')G°(k',iw„). 

Moreover, in the Born approximation, which is appro¬ 
priate away from impurity bound state energies, the mo¬ 


mentum dependence of the potential separates from that 
of the Greens functions convolution when calculating the 
Fourier transformed density of state^^. 


B. Local density of states and impurity types 

Here we will describe how we calculate the various 
LDOS patterns in the work to come. We are consid¬ 
ering FT-STS experiments where the tip of the STM is 
normal (not magnetic or superconducting) and measures 
the LDOS, theoretically this quantity is given as follows 

n{uj, r) = - 7 :^Im (Tr [(1 -1- r^)G(r, r,u; + ir])]) (8) 

ZTT 

where acts on particle-hole degrees of freedom, r] is 
an infinitesimal and G(r, r', cj + ir]) is the double Fourier 
transform of G(k, k',C(; + irj) defined in Eq. 0 above. 
Fourier transforming one finds that the modulation to 
the local density of states in the Born approximation is 


Sn{u,q) = -^;^^Im(Tr [(l-l-T^)G^(k,cj + i7/)E^^^G^(k + q,cj + i7?)]) (9) 


k 


where lm{f{w + ir])) = {f{w + irf) — f{w — ir]))j2i. 

Following Ref. [20] , we are interested in a more general 
interference pattern. Here we imagine that the STM tip 
can resolve different degrees of freedom {e.g. spin can be 
resolved by a spin polarized tip). It is therefore possible 
to insert a general matrix Va into the trace. This matrix 
acts to resolve the component of the LDOS pattern of 
interesfP^. 

We now discuss the various types of perturbations that 
will be considered. Here we will focus on a physical de¬ 


k 


The matrices are outlined in the appendix. The first 
label, a, denotes the type of STM tip. a = 0 is a normal 
(charge) tip while a = 1,2 or 3 resolves the component 
of the electron’s spin along the x^y or z direction respec¬ 
tively. Meanwhile, /3 labels the type of impurity we are 
considering, = 0 is a charge impurity while /d = 1, 2 or 
3 refers to a magnetic impurity with its spin along the x, y 
or z axis. When we refer to patterns such as 
it is the patterns above to which we are referring. It 


scription of these impurities and leave the formal details 
to the Appendix. We begin with a simple charge impu¬ 
rity. Here we imagine that the chemical potential at a 
certain site has been altered. Second, we are interested 
in magnetic impurities. These alter the local Zeeman 
splitting on a single site and so couple to the spin of the 
electrons. In what follows we will combine charge and 
magnetic impurities into a single heading which we will 
refer to as impurities. It is then useful to define the QPI 
patterns for any of these impurities as 


■^)G°(k,w + ir?)F^G°(k + q,w + i7?)]) . (10) 


should be noted that any physical QPI pattern due to 
charge/spin scattering can be written as a combination 
of the above. 

The second class of perturbation we are interested in 
is a local suppression in the superconducting order pa¬ 
rameter. We will refer to this as OP suppression. This 
perturbation can be thought of as a simplified description 
of a vortex where only the OP suppression at the vortex 
core is taken into accounfP^. Alternatively, other types of 
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disorder are sometimes accompanied by OP variationP^^ 
We refer to the LDOS in this case as 6nop- 


C. Model Hamiltonians 

To discuss the surface of a 3DSTI we will make use of 
two model Hamiltonians. The first is a continuum Dirac 
cone model which will be our primary focus. This model 
is a single Dirac con^^ with a chemical potential and 
5-wave pairing. Symbolically we have 


Ho = Hti + Hsc^H^ (11) 

with 

Hti = v j (flic^y.^k^ay - kyax)cy: (12) 

a. 

Hsc = Ao / (ct-|-cLk 4 + C-k^Ck.t) 


where hk = 2X{smkxO-y — sinkyax) — /icrg and Mk = 
e —2t(cos kx^cos ky) and i?k = | (e + 2t(cos k^ + cos ky)) 
where e = 4t. The matrix above is an 8 x 8 matrix, the 
sub matrices hk, Mk and R\^ act on a 2 x 2 space of spin. 
The basis of the above is [(t, 1, T), (|, 1, T), (t, 2, T), (| 
, 2, T), (t, 1, B), (4,, 1, B), (t, 2, B), (i, 2, B)] where arrows 
refer to spin, the numbers refer to one of two bands and 
’T’ and ’H’ refer to top and bottom edges. With this 
basis in mind, we see that Mk couples orbitals within 
the same surface and i?k couples different surfaces. 

We now introduce superconductivity pairing only 
within the same orbital and the same surface. 


A(k) = 


/ Ak 
0 
0 

V 0 


0 

Ak 

0 

0 


0 

0 

Ak 

0 


0 \ 
0 
0 

Ak/ 


(14) 


where Ak = (icr^)Ak. We assume s-wave pairing and 
therefore take Ak = Aq. The above has 8 doubly degen¬ 
erate eigenvalues: 


where Ck = (ck,t 5 Ck 4 )^, /i is the chemical potential 
and Ao is the OP amplitude. We have adopted units 
such that h = 1. Such a model can be diagonal¬ 
ized and yie lds the two eigenva lues at each wave vector 
^±(k) = ^Aq + (±'^|k| — /i)^. Further, finding an ex¬ 
act, closed-form expression for the Green’s function of the 
above system is possible. We have outlined the details 
this calculation in the Appendix. 

The second model is a lattice model of a 3DSTI where 
only the surfaces of this system are considered (we will 
arbitrarily refer to the two surfaces as ’top’ and ’bot¬ 
tom’). This model will be used sparingly and will mostly 
be employed as a physical consistency check for features 
we find in the continuum model. 

In order to model the single Dirac cone on the surface 
of a 3DSTI we follow Marchand and FranJ^. Their ap¬ 
proach is briefly outlined here. The doubling theorem 
states that in a time reversal invariant periodic systems 
Dirac points appear in pairs. Therefore, any two dimen¬ 
sional lattice model can not have an odd number of Dirac 
points. The 3DSTI is indeed a periodic system. However, 
it avoids the doubling theorem by placing half of its Dirac 
points on each surface. To mimic this we employ a two- 
surface model. This is a minimal way to model a surface 
with an odd number of Dirac points. 

In this approach, one begins with a three dimensional 
model and integrates out the bulk, leaving only the two 
edges of the material. Adopting this model gives our 
clean system’s Hamiltonian: 


h{k) 


/K Mk 0 0 \ 

Mk — hk 2i?k 0 
0 hk Mk 

Vo 0 Mk -hj 


(13) 


E^, = ^A2 + (Bk(si,S2)-M)2 (15) 

wherfPSl Bk(si,S2) = siy^Ck + (\/M^ + -^k + S2-Rk)^ 
with si,S 2 = ±1 and = 4A^(sin^/Ca; + ky). The 
clean Green’s function of this system is not analytically 
tractable and so we rely on numerics for its calculation. 


D. Numerical method 


Here we briefly outline the numerical methods we 
use. We numerically compute the LDOS patterns for 
particle-like impurities as described in Eq. ( p!Q| ). The 
idea is to compute the Green’s function in the absence of 
any impurity from the Hamiltonian using G^(k, icj^) = 
{iujm — ^k)~^- Then, we produce the convolution with 
the impurity V^. Once this is obtained, multiplying it 
with then taking the trace produces the wanted re¬ 
sults. 

Each term in the perturbed Green’s function is there¬ 
fore a convolution of two functions of momentum, with 
each of these functions being an entry in the clean sys¬ 
tem’s Green’s function matrix. In order to minimize run 
time and obtain high resolution LDOS maps we preform 
the convolution using the fast Eourier transform algo¬ 
rithm (EET). This amounts to first using EET to express 
the two functions in real space, performing a direct prod¬ 
uct and then using EET again to get back to momentum 
space. This lowers the run time from o(n^) to o(n log n) 
where n is the number of points in the Brillouin zone. 
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III. QUASIPARTICLE INTERFERENCE 
PATTERNS IN A TOPOLOGICAL 
SUPERCONDUCTOR 


A. Zero Chemical Potential 

At the particle-hole symmetric point, i.e., /i = 0, the 
spectrum possesses Lorentz invariance. This enables us 
to find a closed form solution to the continuum model, 
for all of the perturbations we consider. This is done us¬ 
ing the exact Green’s function calculated in the appendix 
along with a standard Feynman parameterizations trick 
(see e.g [28]) which we will not repeat here. For the im¬ 


purities these patterns are most easily presented by defin- 
ing the matrix ^o) = 

the variables zi = 2iujmlQ and Z 2 = 2Ao/^ and ^ = 
\/—zf -h Z 2 and the functions 


F{z) 

G{z) 


2^ 


— 1 — z‘^ arctan 



arctan 

1 

. \/—1—. 

V- 

1 — z‘^ 


(16) 

(17) 


and we have set v = 1 hereafter for simplicity. With 
these definitions we find 


Aa,p{iLUm, Ci,Ao) 



V 



0 

0 

0 


0 

qU2Fz^G{z))-l 
qxQy (1 + z‘^G{z)) 
iqxZiG{z) 


0 

qxqy (1 + z^G(z)) 
ql{2^z^G{z))-l 

iqyZiG{z) 


0 





-iqxZiG{z) 

—iqyZiG{z) 


zl)G{z) 


where Q is an ultraviolet cut-off and qi = qi/\q\. We 
note that with some small exceptions the above results 
are almost identical to those for a strong topological in¬ 
sulator without superconductivit}Eniwith the replacement 
(icOm)^ -> - A^. 

The OP suppression can also be calculated in closed 
form at /i = 0 and yields the results 


Snop{iuJm,ci) 


2Ao{iuJm)G{z) 




(19) 


|k| = — Aq. Our results for the LDOS therefore 

show maximum response at twice this wave vector length. 
Physically, this corresponds to scattering across the di¬ 
ameter of the contours of constant ener gy a s is expected 
for such a circularly symmetric problenP^. This change 
in the singular value of |q| with the frequency uj can be 
seen in the bottom plots of Fig. Here we plot the 
LDOS for several perturbations and probes. The singu¬ 
larities (or kink in the case of Snu) clearly vary with 

UJ. 


For simplicity we have presented to above results in Mat- 
subara space, the proper analytic continuation must be 
employed to obtain the physical QPI patterns. 

We have compared the exact results above with those 
we have obtained from the numerical calculation at /i = 
0 using the model in Eq. (13) with superconductivity 
added. We find that the two patterns agree remark¬ 
ably well showing similar angular and radial features (the 
specifics of these features will be discussed below). We 
have showcased five of these patterns in density plots 
in the top of Fig. We now move on to discuss three 
properties of the QPI patterns we have found in the con¬ 
tinuum model, keeping in mind the agreement between 
the lattice and continuum results. 

First, let us discuss the radial features of the QPI pat¬ 
terns. Inspecting the above results we notice that the 
function G{z) is singular at a q length of g = 2^/uF^^A^ 
while F{z) has a kink at this value (its derivative should 
be singular). To understand this value let us recall 
that the bulk energy b ands for the con tinuum model 
are given by Ek,± = y^Aq + (±|k| — /i)^. With /i = 0 
the two bands become degenerate and contours of con¬ 
stant energy for a given frequency are then given by 


Next, we to study the effect of superconductivity on the 
QPI patterns. In short, the occurrence of Aq in the crit¬ 
ical radius of the LDOS described above is a signature of 
superconductivity. Without superconductivity the peak 
would occur at 2uj and would thus persist all the way 
down to cj = (P. With superconductivity present all fea¬ 
tures disappear once uj < Aq^ i.e. when we probe energy 
scales within the gap. Thus having the pairing present 
in the system shifts the radius of this major feature to 
lower values. We have explored this result in Fig. iby 
plotting several different LDOS patterns for varying uj. 
For u; > Aq we see that the singularity/kink in these pat¬ 
terns moves to smaller values of |q| as uj decreases. For 
cj < Aq, or the solid curve in this figure, we see no signal 
at all. 


Finally, we discuss the angular features of the QPI pat¬ 
terns we have calculated. Studying Eq. (18), we see that 
one requires the input a to be non-zero to find a pat¬ 
tern that is not circularly symmetric. In the case a = 0 
we have either zero (for /3 = 1,2,3) or a circularly sym¬ 
metric function function when = 0. Recall, the (5noo 
result represents a non-magnetic impurity probed with a 
normal STM tip. The other a = 0 results above would 
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FIG. 2: Plots of the LDOS for both the lattice and continuum model with // = 0. Along the top we have subhgures (a) through 
(e) for SriQo, Snu, Sni 2 , Snss, and 5nop (respectively) for the lattice problem. These calculations were done on a 600 x 600 
lattice with A = l,t = 0.5, e = 2, Ao = 0.4 and uj = 0.6. On the bottom subhgures (f) through (j) show plots of the expressions 
for 6noo, Srin, Sni 2 , Sn 33 , and Snop (respectively) in the Dirac model. We have also explored the dependence of these patterns 
on the frequency u. We have hxed Ao = 0.4 and v = 1 and plotted u = 0.385 as a solid line (black online), u = 0.405 as a 
dotted line (red online), uj = 0.425 as a dashed line (blue online) and u = 0.445 as a dashed-dot line (teal online). With the 
exception of Sni 2 , which is plotted along the line qx — qy^ ah plots are along the line qy — 0. 


be those that are obtained from other perturbations and 
a normal STM tip. Thus our underlying chiral system 
does not show any angular dependence along the circu¬ 
lar singularity in the QPI maps. This is the case for a 
normal STM tip at /i = 0 regardless of the perturbation 
type. That being said, as in Ref. m we can see angular 
dependence in the a ^ 0 patterns, which corresponds to 
a spin-filtered STM tip. To look at these angular fea¬ 
tures we have presented density plots of our calculations 
on the lattice model in the top of Fig. 


B. General Chemical Potential 

Moving to non-zero chemical potential makes calcu¬ 
lating the LDOS exactly for the continuum model in¬ 
tractable. Therefore we must rely on numerics for both 
the continuum and lattice model. Our results for the 
/i 7 ^ 0 lattice model are presented in Fig. These re¬ 
sults lead us to make three observations. First, in the 
figure we notice that in addition to one major circular 
pattern, we can see a very subtle secondary pattern in 
the case of magnetic impurities. Second, all patterns es¬ 
sentially respect the same angular symmetry as at /i = 0. 
Third, the radius of the major pattern in the impurity 
QPI patterns changes noticeably when compared to the 
/i = 0 results whereas the OP suppression pattern does 
not. 

Let us put these observations onto some more solid 
ground. We begin with the observation that the angu¬ 


lar symmetry does not change. We have thus established 
that tuning the chemical potential away from the Dirac 
point does not lead to angular signatures in either the 
impurity or OP suppression perturbation. We now ded¬ 
icate the rest of this subsection to explaining the other 
two observations above. 

To understand the above /i dependence, let us recall 
an observation from the last subsection; the major ra¬ 
dial features in the LDOS patterns appear at values of 
q corresponding to impurity scattering between states at 
the same quasiparticle energy. Additionally, this scat¬ 
tering occurs across the diameter of the contour of con¬ 
stant energy. At /i = 0 we have degenerate contours of 
constant energy. When /i is nonzero both of our models 
develop multiple contours of constant energy. In the con¬ 
tinuum mod el, the ra dii of these contours are given by 
|k| = l/i ± ~ ^ol- With these two contours one can 

imagine 4 quasiparticle scattering processes across the 
diameter of these circles. These processes are illustrated 
in Fig. and labelled Ki through K 4 . 

With these different processes in mind we make the fol¬ 
lowing two observations: (1) magnetic and non-magnetic 
impurities favor scattering processes between the same 
contour (intra-contour), (2) the OP suppression pertur¬ 
bation favors scattering processes between different con¬ 
tours (inter-contour). The QPI patterns of impurity 
scaterrers show a major singularity at |q| = 2 Ki for 
/i > 0 with a minor singularity at |q| = 2 K 2 . For /i < 0 
the major singularity occurs at |q| = 2 K 2 with the sec¬ 
ondary minor singularity at |q| = 2 Ki. Meanwhile, the 
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FIG. 3: Plots of the LDOS for the lattice model with /x 7 ^ 0. Subfigures (a) through (e) label (respectively) results for (5noo, 
(5nii, Sni 2 , Sn 33 , and Snop for the lattice problem . These calculations were done on a 600 x 600 lattice with A = 1 , t = 0.5, 
e = 2, Ao = 0.4 and uj = 0.6, fi 7 ^ 0.5 



FIG. 4: Various quasiparticle scattering processes which may 
contribute to the radius of the singularity in the LDOS. The 
broken lines represent the two contours of constant energy 
and the arrow shows the possible process. 



pattern of a n OP sup pression shows singul arities at 2K^ 
when fi < and 27^4 when /jl > “ ^0 

In relation to the above we find interesting math¬ 
ematical relations for the radii of the singularities in 
the LDOS patterns. The impurity QPI patterns have 
their major singularity at |q| = 2 /i + — Aq = 

while the OP supp ression pattern shows a singularity at 
|q| = 2y/c(;^ — Aq = Thus the intra-contour process 
increases linearly with fi while the inter-contour scatter¬ 
ing wave vector length is independent of ja. 

We can understand the above results by heuristic ar¬ 
guments. The two constant energy contours seen in the 
/i 7 ^ 0 case are the result of band reflection as depicted 
in Fig. In this figure the Dirac cone is intermitted 
by a gap and reflected about the e/c = /i horizontal 
line due to the introduction of superconductivity. On 
each of these contours one can define a chirality angle 
(j)k = arg(/ca, iky). This chirality angle affects the wave 
function in two ways. First, even without superconduc¬ 
tivity, spin-orbit coupling locks the spin direction to the 
momentum direction. Second, the pairing inherits this 


FIG. 5: Dispersion with chirality indicated by arrows both be¬ 
fore (left) and after (right) superconductivity is considered. 


chirality angle, becoming effectively a p-wave supercon¬ 
ductor. Its order parameter winding appears in the co¬ 
herence factors. The notion of chirality is depicted by 
arrows in Fig. 

To describe this analytically we note that in our choice 
of basis the Hamiltonian has the form 

( —/i —vke'^^'^ 0 Aq \ 

^ _ —vke~'^^^ —fi —Aq 0 

0 -Ao li -vke-^^^ 

y Aq 0 —vke^^^ ji j 

( 20 ) 

and the positive energy eigenvectors have the form 


IV’±(k)) = 


Te 


^/2 

i 0 k + 


x/2 


* 0 k 

1 

1 


■E’fc,± = 




( 21 ) 



































here Uk,± = = 

—fi ± 'z;|k|. 

The two branches of energy above (denoted with a + 
or a —) have chirality directions which wind in opposite 
directions in the Brillouin zone (see Fig. |^. Notably, 
the chirality direction for a given branch is completely 
inverted upon sending k ^ — k. This has important 
implications for our system where (as we have already 
discussed) scattering across the diameter of contours of 
constant energy is favored. 

Looking at the impurity potentials and considering its 
action on |' 0 ±(k)) one can see that it reverses the chirality 
angle. This means that V'^\'ip±(k)) will have the same 
chirality angle as |'0±(—k)). This explains the favoring 
of the Ki and K 2 transitions. 

Meanwhile, the OP suppression potential replaces both 
the spin and particle hole degree of freedom. As a result, 
the wavefunction is almost in tact (except for the Bogoli- 
ubov coherence factors being exchanged). This means 
V^^\ 2 p±(k)) will have the same chirality angle as |' 0 ±(k)) 
which leads to a suppression of scattering across the same 
contour as here the initial state and the final state have 
opposite chirality angles. Furthermore, these results lead 
to enhancement of inter-contour scattering. This is be¬ 
cause there is a place on the other contour with the same 
chirality angle. For energies above the Dirac point we are 
cutting each branch only once and so the two contours 
we are left with have chirality angles winding in opposite 
directions. This leads to scattering processes. For 
energies below the Dirac point we cut the same energy 
branch twice and so we will have two contours with chi¬ 
rality angles winding in the same direction. For this case 
K 4 processes are favored. 

We explore the above results in our continuum model 
in the left of Fig.|^ Here we have plotted 6 nss and 6 nop 
for several different values of the chemical potential. We 
see very clearly that the impurity scatterers {e.g. the 
Sriss pattern) have singularities that occur at |q| values 
that increase linearly with ji. Meanwhile, the pattern for 
the OP suppression in this figure illustrates clearly that 
the singularity in this pattern does not depend on fi. 

We now cont end that t he functional form with respect 
to fi oi qa = — Aq, the radius for the OP suppres¬ 
sion, and qp = — Aq, the radius for the regular 

impurity, (z.e. linear and constant respectively) is a re¬ 
sult of the underlying Dirac band structur e. To argue this 
we approximate the spectrum hy E = y^Aq + (ek — /x)^ 
where ek is the dispersion of the underlying band struc¬ 
ture. If we wish to invert this equation to find k as a 
function of E, i.e. to find the equation f or the con tours 
of constant energy, we obtain ek = ± \/E‘^ “ ^o* 

Now let us think about a circularly symmetric disper¬ 
sion for simplicity. Then ek = e/^. Let us assume that 
we are close enough to the center of the Brillouin zone 
so that the leading order term in an expansion of ek is 
valid. We then take ek = ek^ where 7 is some number 
and e is a constant. Then the circles of constant energy 


\ I- — 

/ 1/ i/ A 

•j 'll 1 / / 1 / 


, 1 , , , , 1 i , , , i , , , i 1 , , i 


D 0.5 1 1.5 2 

Qx 

2.5 3 

; 


^ J 


‘**'^*.*N 

N \ 

f 

L V 


" \ 

' I ... I ... I ... I ' i 

. I . . . 



0.6 

Q.X 



0.1 

0.3 

0.5 

0.7 

0.9 


FIG. 6: Local density of states for finite chemical potential. 
On the top we show along the qx axis for the continuum 
model (left) and the lattice model (right), while on the bottom 
we show 6n for an anomalous impurity along the qx axis for 
the continuum model (left) and the lattice model (right). In 
both plots we have fixed uj — 0.6, v — Aq = .4 and the 
calculations were done on a 1500 x 1500 lattice of points. In 
each figure we have included vertical dashed lines at relevant 
qx values. In the top figures we have one at 2/i + 2^JZEE^E^ 
for each p considered while at the bottom we have a single 
line at 2\J{jp — Aq. On the bottom we have included a legend 
that labels the value of p used to calculate the data in each 
curve. 


have a radius given by 

We then see that, provided y/E‘^ — Aq and p are of com¬ 
parable size, the only type of dispersion that yields en¬ 
ergy contours that depend on /i in a linear matter is 
7 = 1, or a Dirac-like dispersion. Thus this linear scaling 
of the major singularity in the QPI patterns is a property 
of a linearly dispersing band structure. Furthermore, the 
independence on p of scattering from one contour to an¬ 
other is also a consequence of having a linear dispersion. 
The radius of singularities for such a process will be 

(23) 

The above is p independent for 7 = 1 only. 

The above argument of course only holds provided the 
k values we are interested in are suitably close to the 
F-point (or, more generally, wherever the Dirac point oc¬ 
curs) so that our linear approximation is valid. That is 
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Qx Qx 

FIG. 7: Local density of states for a system in a capacitor. In 
each figure the solid line is the QPI pattern on the top of the 
sample while the broken line is the pattern on the bottom of 
the sample. Along the top we have plotted Snop (left) and 
Snss (right) for V = 0.1 while along the bottom we have Shop 
(left) and Snss (right) for V = 0.4. In both plots we have fixed 
uj = 0.6, r' = 1, Ao = .4, = 0 and the calculations were done 

on a 600 x 600 lattice of points. The plots show cuts of the 
QPI pattern along the line Qy = 0. The vertical dashed line 
in each hgure is a guide to show where the anomalous pattern 
is peaked. 


In our lattice model the two surfaces (top and bottom) 
are only weakly coupled and so we may think of this as in¬ 
troducing a chemical potential to the two surfaces. This 
chemical potential is equal in magnitude but opposite in 
sign on each edge. Let us, without loss of generality, as¬ 
sume the positive bias is introduced on top of the sample 
(notice that Hcap enters the Hamiltonian different by 
a minus sign from the chemical potential). Thinking of 
the two edges as isolated Dirac points, our findings thus 
far in the paper then simply predict the following: The 
radius of the major peak in the LDOS pattern for an im¬ 
purity QPI should decrease on the top of the sample (as 
if /i = —H), and decrease on the bottom (as if /i = H). 
Meanwhile, an OP suppression should look the same on 
the top and bottom of the sample. 

In addition to the qualitative prediction made above, 
the strength of the bias in this setup is, at least in prin¬ 
ciple, tuneable via external means. 

In Fig we present the LDOS on the top and bottom 
of a sample for several values of this bias voltage. We 
see that the exact scenario described above plays out, 
namely as V is increased the singularity in the impurity 
QPI pattern moves to larger |q| on one edge and smaller 
|q| on the other. In addition, we see that the OP sup¬ 
pression placed on the top and bottom of the sample give 
identical results and are essentially independent of this 
bias voltage. 


to say, in a realistic system Ck will contain other sublead¬ 
ing contributions on top of the linear Dirac like term. To 
explore how well our results hold in such a system we 
have calculated 6 nss and 6 nop for several different val¬ 
ues of the chemical potential in our lattice model. These 
results are presented on the right of Fig. In this fig¬ 
ure we see that our observation holds very well provided 
that fi is kept small enough. For intermediate values of 
/i deviations increase as k moves away from the Dirac 
point. 


C. System in a Capacitor 


In the previous subsection we have outlined an inter¬ 
esting dependence of the radii of the major feature in the 
LDOS on the chemical potential. The chemical poten¬ 
tial can be tuned in STIs by chemical doping or, in a 
thin flake, by electrostatic gating. To illustrate the effect 
of a non-zero chemical potential in a simple setting we 
imagine placing a system inside a capacitor which has an 
effect of biasing the system so that there is a potential 
V at the top and — H at the bottom. As mentioned in 
the introduction of this work such a scenario could be 
realized through the use of dual-gate structureP^. We 
describe this by adding 


^Cap = Hdiag(cro, ctq, -ctq, -ctq, -ctq, -ctq, ctq, ctq) (24) 


to our clean Hamiltonian in Eq. (13). 


IV. CONCLUSIONS 

We have studied the QPI patterns induced by differ¬ 
ent local perturbations in both a lattice and a continuum 
model for a SC surface of a three dimensional topologi¬ 
cal insulator. Our results for a half filled band (such that 
the chemical potential is at the Dirac point) are similar 
to those calculated in the past for a strong topological 
insulator surface. The existence of superconductivity in 
the system gaps out the low-lying excitations and shows 
up qualitatively in the radius of the singularity/kink that 
occurs in the QPI pattern. The presence of Aq reduces 
the critical radius of this singularity. For STM bias val¬ 
ues below the superconductive gap no signal can be ob¬ 
served. For non-magnetic impurities the QPI response is 
weak and consists of an edge similar to the normal state 
result. Remarkably though, we find that when disorder 
in the order parameter amplitude is included (one that 
will generically be present) the edge transforms into a 
peak, which should be more easily observed in experi¬ 
ment. The shape of the singularity - edge vs. peak - can 
thus be used as an indicator of the dominant source of 
quasiparticle scattering in the sample. 

With a finite chemical potential we find almost no 
change in the angular features of the QPI pattern. We do 
however find that the quasiparticle scattering processes 
contributing to impurity and OP suppression scattering 
are different. As a result, the singular features of the 
impurity QPI pattern depend linearly on the chemical 
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potential and that those of the OP suppression are inde¬ 
pendent of the chemical potential. We argued that this 
functional dependence is unique to a Dirac-like spectrum 
and showed that it approximately holds in our more com¬ 
plicated lattice model. We have also proposed and veri¬ 
fied with a calculation an alternative method to tune the 
chemical potential by placing the sample in a capacitor. 
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Appendix A: Impurity Potentials 

Here we will describe the mathematical expressions we 
have used to describe the impurity potential. Starting 
with the charge impurity we define _ x<^Tz 

where X = (g) X acts on spin and any other degree 

of freedom in the system and r acts on Nambu space. 
The operator X acts, in a suitable manner, on any other 
degrees of freedom that may be present {i.e. not spin or 
particle-hole). For example, in our lattice model X^^^a' 
acts on the top-bottom surface degree of freedom while 
in the continuum model there are no other degrees of 
freedom left. In the lattice model we consider an impurity 
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localized on a particular edge and so is diagonal 

with entries of 1 on the impurity edge and 0 on the other 
edge. 

Secondly, we are interested in magnetic impurities. 
These alter the local Zeeman splitting on a single site 
and so couple to the spin of the electrons. We con¬ 
sider three separate impurities (one for each cartesian 
direction) while noting that a general magnetic impu¬ 
rity can be written as a linear combination of these po¬ 
tentials. As such we define (in Nambu space) = 

diag where 

Finally, we consider our treatment of the OP suppres¬ 
sion perturbation. This is modelled by taking = 

(icr^ (8)X) (g) (ir^) where we remind the reader that a* act 
on spin degrees of freedom, r* on Nambu space and X is 
as defined above. We calculate the OP QPI pattern by 
using this potential in Eq. 

Let us close this subsection with a brief description of 
the role of the probe potential, W- In the main text we 
have imagined an experimental set-up where the physical 
STM tip is capable of resolving the component of the elec¬ 
tron spin along a particular projection. In this case we 
are interested in only certain components of the change 
in the Greens function matrix. To find these relevant 
components and how they contribute to an interference 
pattern we must place an operator in the trace in Eq. § 
which acts to find the proper contribution. Eor exam¬ 
ple, if we are interested in an STM tip which resolves the 
^-component of the electron spin we would add a to 
Eq. §. In general, we call this matrix where is 

defined in the same way as above. Eor a = 0 we are 
considering only a normal STM tip, while for a = 1, 2, 3 
we consider an STM tip capable of resolving the spin in 
the X, y, or ^ direction respectively. In either case we 
additionally consider the physically relevant case of the 
tip only resolving degrees of freedom on one edge of the 
system; the tip can only make contact with one edge or 
the other. 


Appendix B: Exact Green’s Function for Continuum 
Model 

Here we will outline our calculation of the exact, clean 
Green’s function for the simple Dirac continuum model. 
To find the Green function of Hq we begin by diagonal¬ 
izing Hti- This is accomplished by making the transfor¬ 
mation Ck = where 

_ 1 f _ 

= 1 ) 

where (/)k is the phase of ky ikx • In the above the 
operators 6k,t (^k,t ) annihilate electrons from the upper 
(lower) Dirac cone. Using this transformation we can 
rewrite T/^k = TkT^k where 77 k = (6k, ^Lk)^ 

Tu = ('J- g-? J (B2) 


In terms of this transformation the Green’s function can 
be written 


G(k, iujm) = TkG(k, iWm)Tl (B3) 


where G(k, icj^) is the Green function in the 77 k basis. 
After making this change of basis we can write Hq as 


Ho = 

k 


(B4) 


where 

/ ek.+ 0 -Aoe-*^'< 0 \ 

y ^ 0 ek,- 0 Aoe-*^>' 

-Aoe*‘^>‘ 0 -ek.+ 0 

V 0 Aoe-*^'^ o’ -ek.- / 

(B5) 

where ek,± = —jiXvk. We see from the above matrix 
that the two bands are completely decoupled from each 
other and we have p-wave pairing on each Dirac cone. 
We have the two independent systems 


^k,+ 

Tik,- 


f ek,+ -Aoe 
V-Aoe'^'' -ek,+ ) 

( ek.- Aoe-*^'<\ 
l^Aoe*^>' -ek.- ) 


(B6) 


Defining the two Green’s function {iuJrn — ^k,A)Gk,A = 1 
where A = ±1 it is straightforward to show 

A 1 (+ ek.A -AAoe“®'^''\ . . 

{ioJrtiY - -Ck.A “ ^k.A ) 

where E\^ ± = + Aq. Our Full Green function is 

then 


G\z.{iuJrn) 

where 


/ 5k(*Wm) -Aoe fu{i(^Tn)\ 

(B8) 


~l~€k,+ 



igjyn, -|-Ck, — 

0 

- -El 


(B9) 


Now we apply the matrices Tk to get back to the Green’s 
function in a spin basis. It reads 


Gk(7'^m) 


/^k(^^ m m) \ 


(BIO) 


where g-k{iu:m) = t^k 5 k(*Wm)[/,[ and fk{ii^m) = 

6/kpk(^<^m)^Tk- Performing the matrix multiplication 
one can show that 


pk('^^m) 

5k,+ +5k.- (flk,-- 5k,+)A 

2 (sk.- - 5k.+) 5k.+ + 5k.- J 


(Bll) 
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and 


where 5k,A = ^ /k,A 

fkiiiOm) = (B12) 

1/e*‘^'‘(/k,+-/k,-) (/k,-+/k,+) \ 

2 - (/k,- + /k.+) -e-‘^M/k,+ - /k .-)J 


_ ^0 _ 



